package scu.maqiang.elastodynamics;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L2PE;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class PEDynamics1DL2 {

	public static void main(String[] args) {
		Mesh1L2 mesh = new Mesh1L2();
		mesh.line(0, 1, 128);
		
		FES1L2PE fs = new FES1L2PE(mesh);
		//FES1L21 fs1 = new FES1L21(mesh);
		//SRMatrix MM = new SRMatrix(fs1.GetNdof());
		SRMatrix K = new SRMatrix(fs.GetNdof());
		SRMatrix M0 = new SRMatrix(fs.GetNdof());
		SRMatrix M = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		double[] u = new double[fs.GetNdof()];
		double[] v = u.clone();
		double[] a = u.clone();
		double[] u0 = u.clone();
		double[] v0 = u.clone();
		double[] a0 = u.clone();
		
		double time = 0.0, dt = 0.005, endTime = 10.0;
		double alpha = 0.25, delta = 0.5;
		double c0 = 1.0 / (alpha * dt * dt);
		double c1 = delta / (alpha * dt);
		double c2 = 1.00 / (alpha * dt);
		double c3 = 0.50 / alpha - 1.0;
		double c4 = delta / alpha - 1.0;
		double c5 = 0.5 * dt * (delta / alpha - 2.0);
		double c6 = dt * (1 - delta);
		double c7 = delta * dt;
		
		double G1 = 1.0, ep1 = 0.5, e1 = 0.1;
		fs.assemblePEStiff(new double[] {G1, ep1, e1}, BVPType.COMMON, K);
		fs.assemblePEMass(new double[] {1.0}, BVPType.CONSIST_MASS, M0);
		fs.assemblePEMass(new double[] {1.0}, BVPType.CONSIST_MASS, M);
		//fs1.assembleHeatMass(new double[] {1.0}, BVPType.CONSIST_MATRIX, MM);
		//System.out.println(M);
		SRMatrix L = K.add(c0, M);
		double[] uva = u.clone();
		//���ó�ʼλ��
		ScalarFunc initFunc =  (x, label, param) -> Math.sin(0.5 * Math.PI * x[0]); 
		double[] uu0 = fs.valueFromFunc(initFunc, null);
		
		for(int i = 0; i < uu0.length; i++) {
			u0[2 * i] = uu0[i];
		}
		double[] Ku0 = K.mul(u0);
		//K.sort();
		//System.out.println(K);
		//DirectSSolver mSolver = new DirectSSolver(M0, fs.getnBand());
		//M0.sort();
		//System.out.println(M0);
		//System.exit(0);
		for(int i = 0; i < fs.GetNdof(); i++) {
			RHS[i] = RHS[i] - Ku0[i];
		}
		
		//mSolver.PLUSolve(RHS, a0, 1);
		Arrays.setAll(a, i -> a0[i]);
		System.exit(0);
		
		for(int i = 0; i < fs.GetNdof(); i++) {
			uva[i] = c0 * u0[i] + c2 * v0[i] + c3 * a0[i];
		}
		M.mul(uva, RHS);
		String outputFile1 = "PEFreeVibration1Du.dat";
		String outputFile2 = "PEFreeVibration1Dv.dat";
		time = 0.0;
		
		double[][] uvtime = new double[2][mesh.getNv()];
		fs.extractComponent(u0, uvtime);
		mesh.toTecplot(outputFile1, time, uvtime[0]);
		mesh.toTecplot(outputFile2, time, uvtime[1]);
		fs.applyBC_MBN(L, Direct.X, 1);
		fs.applyBC_MBN(L, Direct.Y, 2);
//		L.sort();
//		//System.out.println(L);
//		System.out.println(L.getNNZ());
		//����ʼ���ٶ�
		
		
		NewIterSSolver solver = new NewIterSSolver(L);
//		ScalarCoefFunc BCValueFunc = (xy, llabel, param) -> 0.1 * Math.sin(8 * Math.PI * param[0]);
		
		//myParam[0] = 0.005;
		//fs.applyBC_MBN(RHS, Direct.X, new RegionFunc[] {BCValueFunc}, myParam, exteriorLabel);
		//for(int i = 0; i < fs.GetNdof(); i++) {
		//	System.out.println(RHS[i]);
		//}
//		ScalarCoefFunc[] bdRegionFunc = new ScalarCoefFunc[] {BCValueFunc};

		double[] myParam = new double[1];
		while(time <= endTime) {
			time += dt;
			myParam[0] = time;
			System.out.println("time: " + time);
			//fs.applyBC_MBN(RHS, Direct.All, 0.0, 1, 2);
			fs.applyBC_MBN(RHS, Direct.X, 0.0, 1);
			fs.applyBC_MBN(RHS, Direct.Y, 0.0, 2);
			//System.out.println(Arrays.toString(RHS));
			
			solver.PCGSSOR(RHS, u, 1.5, 1);
			for(int i = 0; i < fs.GetNdof(); i++) {
			    a[i] = c0 * (u[i] - u0[i]) - c2 * v0[i] - c3 * a0[i];
			    v[i] = v0[i] + c6 * a0[i] + c7 * a[i]; 
			    uva[i] = c0 * u[i] + c2 * v[i] + c3 * a[i];				
			}
			M.mul(uva, RHS);
			for(int i = 0; i < fs.GetNdof(); i++) {
				u0[i] = u[i];
			    v0[i] = v[i];
			    a0[i] = a[i];
			}
			fs.extractComponent(u0, uvtime);
			mesh.toTecplot(outputFile1, time, uvtime[0]);
			mesh.toTecplot(outputFile2, time, uvtime[1]);
		}
	}

}
